{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "</center>\n",
    "Автор материала: программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Тема 7. Обучение без учителя: PCA и кластеризация\n",
    "## <center>Бонус. Метод главных компонент. Игрушечный пример"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Пусть дана выборка X.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array([[1.0, 3.0], [3.0, 5.0], [5.0, 1.0], [7.0, 4.0], [4.0, 7.0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Как выбрать направление, в проекции на которое дисперсия координат точек максимальна? Синия прямая или зеленая? А может, красная?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1])\n",
    "plt.plot(np.linspace(1, 8, 10), np.linspace(1, 8, 10))\n",
    "plt.plot(np.linspace(1, 8, 10), np.linspace(2, 4, 10))\n",
    "plt.plot(np.linspace(1, 8, 10), np.linspace(5, 2, 10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Стандартизуем матрицу X. Вычитаем средние по столбцам (4 и 4) и делим на стандартные отклонения по столбцам (2 и 2). Кстати, пришлось писать код, чтоб подобрать координаты так, чтоб все средние и отклонения были целыми :)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled = StandardScaler().fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_scaled[:, 0], X_scaled[:, 1])\n",
    "plt.plot([-2, 2], [0, 0], c=\"black\")\n",
    "plt.plot([0, 0], [-2, 2], c=\"black\")\n",
    "plt.xlim(-2, 2)\n",
    "plt.ylim(-2, 2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Назовем новые координаты (стоблцы матрицы X_scaled) $x_1$ и $x_2$. Задача: найти такую линейную комбинацию $z = \\alpha x_1 + \\beta x_2$, что дисперсия $z$ максимальна. При этом должно выполняться $\\alpha^2 + \\beta^2 = 1.$**\n",
    "**Заметим что $$\\Large D[z] = E[(z - E[z])^2]  = E[z^2] = \\frac{1}{n} \\sum_i^n z_i^2,$$ поскольку $E[z] = \\alpha E[x_1] + \\beta E[x_2] = 0$ (новые координаты центрированы).**\n",
    "\n",
    "**Тогда задача формализуется так:**\n",
    "$$\\Large \\begin{cases} \\max_{\\alpha, \\beta} \\sum_i^n (\\alpha x_{1_i} + \\beta x_{2_i})^2 \\\\  \\alpha^2 + \\beta^2 = 1\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "У нас $2z = [-3\\alpha -\\beta,\\ -\\alpha +\\beta,\\ \\alpha -3\\beta,\\ 3\\alpha,\\ 3\\beta]^T$ (Для задачи максимизации неважно, что мы умножили на 2, зато так удобней).\n",
    "           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Распишем в нашем случае: $  \\sum_i^n (\\alpha x_{1_i} + \\beta x_{2_i})^2 = (-3\\alpha -\\beta)^2 + ( -\\alpha +\\beta)^2 +( \\alpha -3\\beta)^2 +( 3\\alpha)^2 +( 3\\beta)^2 = 20\\alpha^2 - 2\\alpha\\beta + 20\\beta^2$ =  <font color='green'>\\\\ поскольку $\\alpha^2 + \\beta^2 = 1$ \\\\ </font>  = $20 - 2\\alpha\\beta$. Осталось только минимизировать $\\alpha\\beta$. Можно это делать методом Лагранжа, но в данном случае можно проще\n",
    "\n",
    "$$\\Large \\begin{cases} \\min_{\\alpha, \\beta} \\alpha\\beta \\\\  \\alpha^2 + \\beta^2 = 1\\end{cases}$$\n",
    "\n",
    "$\\Large \\alpha\\beta = \\beta^2(\\frac{\\alpha}{\\beta})$ = <font color='green'>\\\\ замена t = $\\frac{\\alpha}{\\beta}, \\alpha^2 + \\beta^2 = 1$ \\\\  </font> = $\\Large \\frac{t}{1+t^2}$. Ищем минимум функции одной переменной, находим, что $t^* = -1$.\n",
    "\n",
    "Значит, $$\\Large \\begin{cases} \\alpha^* = -\\beta^*\\\\  (\\alpha^*)^2 + (\\beta^*)^2 = 1\\end{cases} \\Rightarrow \\alpha^* = \n",
    "\\frac{1}{\\sqrt{2}}, \\beta^* = - \\frac{1}{\\sqrt{2}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Итак, $$\\Large z = \\frac{1}{\\sqrt{2}} x_1 - \\frac{1}{\\sqrt{2}}x_2$$ То есть ось $z$ повернута на 45 градусов относительно $x_1$ и $x_2$ и \"направлена на юго-восток\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_scaled[:, 0], X_scaled[:, 1])\n",
    "plt.plot([-2, 2], [0, 0], c=\"black\")\n",
    "plt.plot([0, 0], [-2, 2], c=\"black\")\n",
    "plt.plot([-2, 2], [2, -2], c=\"red\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Новые координаты точек по оси z:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled.dot(np.array([1.0 / np.sqrt(2), -1.0 / np.sqrt(2)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Сингулярное разложение матрицы X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Представление будет таким: $X = U\\Sigma V^T$. \n",
    "\n",
    "- Матрица $U$ составлена из собственных векторов матрицы $XX^T$. Это левые сингулярные векторы матрицы $X$;\n",
    "- Матрица $V$ составлена из собственных векторов матрицы $X^TX$. Это правые сингулярные векторы матрицы $X$;\n",
    "- Матрица $\\Sigma$ - диагональная (вне главной диагонали нули), и на диагонали стоят корни  из собственных значений матрицы $X^TX$ (или $XX^T$). Это сингулярные числа матрицы $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$XX^T$ выглядит так:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled.dot(X_scaled.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$X^TX$ выглядит так:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled.T.dot(X_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Собственные вектора $XX^T$ (левые сингулярные):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.eig(X_scaled.dot(X_scaled.T))[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Собственные вектора $X^TX$ (правые сингулярные). Эти вектора задают представление главных компонент через исходные координаты (то есть они задают поворот)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.eig(X_scaled.T.dot(X_scaled))[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что главные компоненты: $$\\Large z_1 = \\frac{1}{\\sqrt{2}} x_1 - \\frac{1}{\\sqrt{2}}x_2,\\ z_2 = \\frac{1}{\\sqrt{2}} x_1 + \\frac{1}{\\sqrt{2}}x_2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Собственные значения $X^TX$ (сингулярные числа):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.eig(X_scaled.T.dot(X_scaled))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.eig(X_scaled.dot(X_scaled.T))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import svd\n",
    "\n",
    "U, Sigma, VT = svd(X_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Действительно. На диагонали матрицы $\\Sigma$ стоят корни из собственных значений $X^TX$ ($\\sqrt{5.25} \\approx 2.29, \\sqrt{4.75} \\approx 2.18$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вектора матрицы $VT$ (правые сингулярные векторы для исходной матрицы) задают поворот. То есть первая главная компонента \"смотрит на юго-восток\", вторая - на юго-запад."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Представление данных в проекции на 2 главные компоненты $Z = XV$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled.dot(VT.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_scaled[:, 0], X_scaled[:, 1])\n",
    "plt.plot([-2, 2], [0, 0], c=\"black\")\n",
    "plt.plot([0, 0], [-2, 2], c=\"black\")\n",
    "plt.plot([-2, 2], [2, -2], c=\"red\")\n",
    "plt.plot([-2, 2], [-2, 2], c=\"red\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Здесь SVD SciPy \"направил\" ось z1 вправо и вниз, а ось z2 - влево и вниз. Можно проверить, что представление получилось правильным. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
